##################################################
# Replication Code
# Taeyong Park and Andrew Reeves
# "Local Unemployment and Voting for President: Uncovering Causal Mechanisms"
# Summary: Replicate Tables 1 - 6 in the online appendix.
# RUN MODELS: Code for running all the models in the tables.
# REPORT RESULTS: Code for reporting the results.
# The .bug files must be in the working directory.
##################################################


###########################
#
# RUN MODELS
#
###########################

############################
#
# 2008: Table 1 and Table 4      
#
############################

library(foreign); library(stringr); library(R2WinBUGS);library(rjags);library(R2jags);library(coda);library(superdiag); library(MCMCpack)

rm(list = ls())
data = read.dta("data08Imputed1.dta")

# Evaluations
data$VBad = rep(0, nrow(data))
data$VBad[data$natecon5==5] = 1
data$Bad = rep(0, nrow(data))
data$Bad[data$natecon5==4] = 1
data$Good = rep(0, nrow(data))
data$Good[data$natecon5==2] = 1
data$VGood = rep(0, nrow(data))
data$VGood[data$natecon5==1] = 1

# Age
data$younger30 = rep(0, nrow(data))
data$younger30[data$age<30] = 1
data$older65 = rep(0, nrow(data))
data$older65[data$age>=65] = 1

# Race
data$black = rep(0, nrow(data))
data$black[data$race==1]=1
data$latino = rep(0, nrow(data))
data$latino[data$race==2]=1

# Employment status
data$fullTime = rep(0, nrow(data))
data$fullTime[data$employment == 1] = 1 
data$partTime = rep(0, nrow(data))
data$partTime[data$employment == 2] = 1 
data$unemployed = rep(0, nrow(data))
data$unemployed[data$employment == 3] = 1 

# Education
data$noHighSchool = rep(0, nrow(data))
data$noHighSchool[data$educ==4] = 1
data$HighSchool = rep(0, nrow(data))
data$HighSchool[data$educ==3] = 1
data$someCollege = rep(0, nrow(data))
data$someCollege[data$educ==1] = 1
data$fourCollege = rep(0, nrow(data))
data$fourCollege[data$educ==2] = 1

# Party ID
data$Rep = rep(0, nrow(data))
data$Rep[data$party3==3] = 1
data$Ind = rep(0, nrow(data))
data$Ind[data$party3==2] = 1
data$Dem = rep(0, nrow(data))
data$Dem[data$party3==1] = 1

# News interest
data$Very=ifelse(data$newsInt==4, 1, 0)
data$Some=ifelse(data$newsInt==3, 1 ,0)
data$Few=ifelse(data$newsInt==1 | data$newsInt==2, 1, 0)

# Unemployment
data$avgUnemp07 = apply(data[,7:18], 
                        FUN=mean, MARGIN=1)
data$avgUnemp08 = apply(data[,19:30], 
                        FUN=mean, MARGIN=1)
data$chgUnemp0807=data$avgUnemp08-data$avgUnemp07

# Income 
data$chgInc0807=(data$inc08-data$inc07)/1000
data$inc08=data$inc08/1000

# Gas price
data$avgGas = (data$gasOct08 + data$gasSept08 + data$gasAug08)/3

# Foreclosure
data$totForcAugOct = data$forcAug08  + data$forcSep08+   data$forcOct08


#################
# TABLE 1 
#################

## Model 1 in Table 1 ##
attach(data)
y = natecon5
N = nrow(data)
Nstate = length(unique(data$stateID)) 
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "stateID"=stateID, "chgUnemp0807"= chgUnemp0807, "avgUnemp08"=avgUnemp08, 
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Nstate" = Nstate, "Ncat" = Ncat)

params = c("b.state", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", 
           "b.20", "b.21", "cut", "sigma.state")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2008_M1.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed=0508)
mcmc08_step1Model1_M1 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:21){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc08_step1Model1_M1)[1:21][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc08_step1Model1_M1)[order(mcmcID)] = c("chgUnemp0807", "avgUnemp08",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc08_step1Model1_M1, "mcmc08_step1Model1_M1")
rm(output); rm(mcmc08_step1Model1_M1)


## Model 2 in Table 1 ##
attach(data)
y = natecon5
N = nrow(data)
Nstate = length(unique(data$stateID)) 
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "stateID"=stateID, "chgUnemp0807"= chgUnemp0807, "avgUnemp08" = avgUnemp08, "chgInc0807"=chgInc0807, "inc08"=inc08,
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Nstate" = Nstate, "Ncat" = Ncat)

params = c("b.state", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "b.22", "b.23", "cut", "sigma.state")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2008_M2.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed=0508)
mcmc08_step1Model1_M2 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:23){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc08_step1Model1_M2)[1:23][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc08_step1Model1_M2)[order(mcmcID)] = c("chgUnemp0807", "avgUnemp08","chgInc0807", "inc08",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc08_step1Model1_M2, "mcmc08_step1Model1_M2")
rm(output); rm(mcmc08_step1Model1_M2)


## Model 3 in Table 1 ##
attach(data)
y = natecon5
N = nrow(data)
Nstate = length(unique(data$stateID)) 
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "stateID"=stateID, "chgUnemp0807"= chgUnemp0807, "avgUnemp08" = avgUnemp08, "chgInc0807"=chgInc0807, "inc08"=inc08,
                 "avgGas"=avgGas, "totForcAugOct"=totForcAugOct,
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Nstate" = Nstate, "Ncat" = Ncat)

params = c("b.state", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "b.22", "b.23", "b.24", "b.25", "cut", "sigma.state")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2008.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed=0508)
mcmc08_step1Model1_M3 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:25){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc08_step1Model1_M3)[1:25][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc08_step1Model1_M3)[order(mcmcID)] = c("chgUnemp0807", "avgUnemp08","chgInc0807", "inc08",
                                                   "avgGas", "totForcAugOct",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc08_step1Model1_M3, "mcmc08_step1Model1_M3")
rm(output); rm(mcmc08_step1Model1_M3)


## Model 4 in Table 1 ##
attach(data)
y = natecon5
N = nrow(data)
Ncounty = length(unique(data$countyID)) 
countyID = as.integer(droplevels(factor(countyID)))
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "countyID"=countyID, "chgUnemp0807"= chgUnemp0807, "avgUnemp08" = avgUnemp08, "chgInc0807"=chgInc0807, "inc08"=inc08,
                 "avgGas"=avgGas, "totForcAugOct"=totForcAugOct,
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Ncounty" = Ncounty, "Ncat" = Ncat)

params = c("b.county", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "b.22", "b.23", "b.24", "b.25", "cut", "sigma.county")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2008_M4.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed=0508)
mcmc08_step1Model1_M4 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:25){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc08_step1Model1_M4)[1:25][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc08_step1Model1_M4)[order(mcmcID)] = c("chgUnemp0807", "avgUnemp08","chgInc0807", "inc08",
                                                   "avgGas", "totForcAugOct",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc08_step1Model1_M4, "mcmc08_step1Model1_M4")
rm(output); rm(mcmc08_step1Model1_M4)



#############
# TABLE 4
#############

## Model 1 in Table 4 ##

model08 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp0807 +  avgUnemp08 +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="stateID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc08_step1Model2_M1 = model08$mcmc[, 1:27]
colnames(mcmc08_step1Model2_M1)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp0807", "avgUnemp08",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc08_step1Model2_M1, "mcmc08_step1Model2_M1")
rm(mcmc08_step1Model2_M1)


## Model 2 in Table 4 ##
model08 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp0807 +  avgUnemp08 + chgInc0807 + inc08 +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="stateID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc08_step1Model2_M2 = model08$mcmc[, 1:29]
colnames(mcmc08_step1Model2_M2)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp0807", "avgUnemp08", "chgInc0807", "inc08",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc08_step1Model2_M2, "mcmc08_step1Model2_M2")
rm(mcmc08_step1Model2_M2)

## Model 3 in Table 4 ##
model08 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp0807 +  avgUnemp08 + chgInc0807 + inc08 +
                       avgGas + totForcAugOct +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="stateID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc08_step1Model2_M3 = model08$mcmc[, 1:31]
colnames(mcmc08_step1Model2_M3)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp0807", "avgUnemp08", "chgInc0807", "inc08",
                                  "avgGas", "totForcAugOct",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc08_step1Model2_M3, "mcmc08_step1Model2_M3")
rm(mcmc08_step1Model2_M3)


## Model 4 in Table 4 ##
model08 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp0807 +  avgUnemp08 + chgInc0807 + inc08 +
                       avgGas + totForcAugOct +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="countyID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc08_step1Model2_M4 = model08$mcmc[, 1:31]
colnames(mcmc08_step1Model2_M4)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp0807", "avgUnemp08", "chgInc0807", "inc08",
                                  "avgGas", "totForcAugOct",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc08_step1Model2_M4, "mcmc08_step1Model2_M4")
rm(mcmc08_step1Model2_M4)

## Model 5 in Table 4 ##
data_noVGood=data[data$natecon5 != 1,]

model08 = MCMChlogit(fixed=pvote2 ~ Good + Bad + VBad +
                       chgUnemp0807 +  avgUnemp08 + chgInc0807 + inc08 +
                       avgGas + totForcAugOct +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="stateID",
                     data=data_noVGood, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc08_step1Model2_M5 = model08$mcmc[, 1:30]
colnames(mcmc08_step1Model2_M5)=c("intercept", "Good", "Bad", "VBad", 
                                  "chgUnemp0807", "avgUnemp08", "chgInc0807", "inc08",
                                  "avgGas", "totForcAugOct",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc08_step1Model2_M5, "mcmc08_step1Model2_M5")
rm(mcmc08_step1Model2_M5)


############################
#
# 2012: Table 2 and Table 5      
#
############################

library(stringr)
library(R2WinBUGS);library(rjags);library(R2jags);library(coda);library(superdiag); library(MCMCpack)

rm(list = ls())
data = read.dta("data12Imputed1.dta")

# Evaluations
data$VBad = rep(0, nrow(data))
data$VBad[data$natecon5==5] = 1
data$Bad = rep(0, nrow(data))
data$Bad[data$natecon5==4] = 1
data$Good = rep(0, nrow(data))
data$Good[data$natecon5==2] = 1
data$VGood = rep(0, nrow(data))
data$VGood[data$natecon5==1] = 1

# Age
data$younger30 = rep(0, nrow(data))
data$younger30[data$age<30] = 1
data$older65 = rep(0, nrow(data))
data$older65[data$age>=65] = 1

# Race
data$black = rep(0, nrow(data))
data$black[data$race==1]=1
data$latino = rep(0, nrow(data))
data$latino[data$race==2]=1

# Employment status
data$fullTime = rep(0, nrow(data))
data$fullTime[data$employment == 1] = 1 
data$partTime = rep(0, nrow(data))
data$partTime[data$employment == 2] = 1 
data$unemployed = rep(0, nrow(data))
data$unemployed[data$employment == 3] = 1 

# Education
data$noHighSchool = rep(0, nrow(data))
data$noHighSchool[data$educ==4] = 1
data$HighSchool = rep(0, nrow(data))
data$HighSchool[data$educ==3] = 1
data$someCollege = rep(0, nrow(data))
data$someCollege[data$educ==1] = 1
data$fourCollege = rep(0, nrow(data))
data$fourCollege[data$educ==2] = 1

# Party ID
data$Rep = rep(0, nrow(data))
data$Rep[data$party3==3] = 1
data$Ind = rep(0, nrow(data))
data$Ind[data$party3==2] = 1
data$Dem = rep(0, nrow(data))
data$Dem[data$party3==1] = 1

# News interest
data$Very=ifelse(data$newsInt==4, 1, 0)
data$Some=ifelse(data$newsInt==3, 1 ,0)
data$Few=ifelse(data$newsInt==1 | data$newsInt==2, 1, 0)


# Unemployment
data$avgUnemp11 = apply(data[,15:26], 
                        FUN=mean, MARGIN=1)
data$avgUnemp12 = apply(data[,27:38], 
                        FUN=mean, MARGIN=1)
data$chgUnemp1211=data$avgUnemp12-data$avgUnemp11

# Income
data$chgInc1211=(data$inc12-data$inc11)/1000
data$inc12 = data$inc12/1000

# Gas price
data$avgGas = (data$gasOct + data$gasSep + data$gasAug)/3


#################
# TABLE 2 
#################

## Model 1 in Table 2 ##
attach(data)
y = natecon5
N = nrow(data)
Nstate = length(unique(data$stateID)) 
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "stateID"=stateID, "chgUnemp1211"= chgUnemp1211, "avgUnemp12"=avgUnemp12,
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Nstate" = Nstate, "Ncat" = Ncat)

params = c("b.state", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "cut", "sigma.state")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2012_M1.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed=0508)
mcmc12_step1Model1_M1 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:21){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc12_step1Model1_M1)[1:21][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc12_step1Model1_M1)[order(mcmcID)] = c("chgUnemp1211", "avgUnemp12",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc12_step1Model1_M1, "mcmc12_step1Model1_M1")
rm(output); rm(mcmc12_step1Model1_M1)

## Model 2 in Table 2 ##
attach(data)
y = natecon5
N = nrow(data)
Nstate = length(unique(data$stateID)) 
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "stateID"=stateID, "chgUnemp1211"= chgUnemp1211, "avgUnemp12" = avgUnemp12, "chgInc1211"=chgInc1211, "inc12"=inc12,
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Nstate" = Nstate, "Ncat" = Ncat)

params = c("b.state", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "b.22", "b.23", "cut", "sigma.state")
output = jags(data.jags, inits=NULL, params, model.file="Model1_M2012_M2.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed = 0508)

mcmc12_step1Model1_M2 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:23){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc12_step1Model1_M2)[1:23][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc12_step1Model1_M2)[order(mcmcID)] = c("chgUnemp1211", "avgUnemp12","chgInc1211", "inc12",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc12_step1Model1_M2, "mcmc12_step1Model1_M2")
rm(output); rm(mcmc12_step1Model1_M2)

## Model 3 in Table 2 ##
attach(data)
y = natecon5
N = nrow(data)
Nstate = length(unique(data$stateID)) 
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "stateID"=stateID, "chgUnemp1211"= chgUnemp1211, "avgUnemp12" = avgUnemp12, "chgInc1211"=chgInc1211, "inc12"=inc12,
                 "avgGas"=avgGas, "totForcAugOct"=totForcAugOct,
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Nstate" = Nstate, "Ncat" = Ncat)

params = c("b.state", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "b.22", "b.23", "b.24", "b.25", "cut", "sigma.state")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2012.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed = 0508)

mcmc12_step1Model1_M3 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:25){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc12_step1Model1_M3)[1:25][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc12_step1Model1_M3)[order(mcmcID)] = c("chgUnemp1211", "avgUnemp12","chgInc1211", "inc12",
                                                   "avgGas", "totForcAugOct",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc12_step1Model1_M3, "mcmc12_step1Model1_M3")
rm(output); rm(mcmc12_step1Model1_M3)


## Model 4 in Table 2 ##
attach(data)
y = natecon5
N = nrow(data)
Ncounty = length(unique(data$countyID)) 
countyID = as.integer(droplevels(factor(countyID)))
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "countyID"=countyID, "chgUnemp1211"= chgUnemp1211, "avgUnemp12" = avgUnemp12, "chgInc1211"=chgInc1211, "inc12"=inc12,
                 "avgGas"=avgGas, "totForcAugOct"=totForcAugOct,
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Ncounty" = Ncounty, "Ncat" = Ncat)

params = c("b.county", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "b.22", "b.23", "b.24", "b.25", "cut", "sigma.county")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2012_M4.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed = 0508)

mcmc12_step1Model1_M4 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:25){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc12_step1Model1_M4)[1:25][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc12_step1Model1_M4)[order(mcmcID)] = c("chgUnemp1211", "avgUnemp12","chgInc1211", "inc12",
                                                   "avgGas", "totForcAugOct",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc12_step1Model1_M4, "mcmc12_step1Model1_M4")
rm(output); rm(mcmc12_step1Model1_M4)


#############
# Table 5
#############

## Model 1 in Table 5 ##
model12 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp1211 + avgUnemp12 + 
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="stateID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc12_step1Model2_M1 = model12$mcmc[, 1:27]
colnames(mcmc12_step1Model2_M1)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp1211", "avgUnemp12",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc12_step1Model2_M1, "mcmc12_step1Model2_M1")
rm(mcmc08_step1Model2_5)

## Model 2 in Table 5 ##
model12 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp1211 +  avgUnemp12 + chgInc1211 + inc12 +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="stateID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc12_step1Model2_M2 = model12$mcmc[, 1:29]
colnames(mcmc12_step1Model2_M2)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp1211", "avgUnemp12", "chgInc1211", "inc12",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc12_step1Model2_M2, "mcmc12_step1Model2_M2")
rm(mcmc12_step1Model2)

## Model 3 in Table 5 ##
model12 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp1211 +  avgUnemp12 + chgInc1211 + inc12 +
                       avgGas + totForcAugOct +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="stateID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc12_step1Model2_M3 = model12$mcmc[, 1:31]
colnames(mcmc12_step1Model2_M3)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp1211", "avgUnemp12", "chgInc1211", "inc12",
                                  "avgGas", "totForcAugOct",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc12_step1Model2_M3, "mcmc12_step1Model2_M3")



## Model 4 in Table 5 ##
model12 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp1211 +  avgUnemp12 + chgInc1211 + inc12 +
                       avgGas + totForcAugOct +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="countyID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc12_step1Model2_M4 = model12$mcmc[, 1:31]
colnames(mcmc12_step1Model2_M4)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp1211", "avgUnemp12", "chgInc1211", "inc12",
                                  "avgGas", "totForcAugOct",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc12_step1Model2_M4, "mcmc12_step1Model2_M4")


############################
#
# 2016: Table 3 and Table 6      
#
############################

library(stringr)
library(R2WinBUGS);library(rjags);library(R2jags);library(coda);library(superdiag); library(MCMCpack)
library(foreign)

data = read.dta("data16Imputed1.dta") 

# Evaluations
data$VBad = rep(0, nrow(data))
data$VBad[data$natecon5==5] = 1
data$Bad = rep(0, nrow(data))
data$Bad[data$natecon5==4] = 1
data$Good = rep(0, nrow(data))
data$Good[data$natecon5==2] = 1
data$VGood = rep(0, nrow(data))
data$VGood[data$natecon5==1] = 1

# Age
data$younger30 = rep(0, nrow(data))
data$younger30[data$age<30] = 1
data$older65 = rep(0, nrow(data))
data$older65[data$age>=65] = 1

# Race
data$black = rep(0, nrow(data))
data$black[data$raceNew==1]=1
data$latino = rep(0, nrow(data))
data$latino[data$raceNew==2]=1

# Employment status
data$fullTime = rep(0, nrow(data))
data$fullTime[data$employment == 1] = 1 
data$partTime = rep(0, nrow(data))
data$partTime[data$employment == 2] = 1 
data$unemployed = rep(0, nrow(data))
data$unemployed[data$employment == 3] = 1 

# Education
data$noHighSchool = rep(0, nrow(data))
data$noHighSchool[data$educNew==4] = 1
data$HighSchool = rep(0, nrow(data))
data$HighSchool[data$educNew==3] = 1
data$someCollege = rep(0, nrow(data))
data$someCollege[data$educNew==1] = 1
data$fourCollege = rep(0, nrow(data))
data$fourCollege[data$educNew==2] = 1

# Party ID
data$Rep = rep(0, nrow(data))
data$Rep[data$party3==3] = 1
data$Ind = rep(0, nrow(data))
data$Ind[data$party3==2] = 1
data$Dem = rep(0, nrow(data))
data$Dem[data$party3==1] = 1

# News interest
data$Very=ifelse(data$newsInt==4, 1, 0)
data$Some=ifelse(data$newsInt==3, 1 ,0)
data$Few=ifelse(data$newsInt==1 | data$newsInt==2, 1, 0)

# Unemp
data$avgUnemp15 = apply(data[,5:16], 
                        FUN=mean, MARGIN=1)
data$avgUnemp16 = apply(data[,17:28], 
                        FUN=mean, MARGIN=1)
data$chgUnemp1615=data$avgUnemp16-data$avgUnemp15

# Income
data$chgInc1615=(data$inc16-data$inc15)/1000
data$inc16 = data$inc16/1000


#################
# Table 3 
#################

## Model 1 in Table 3 ##
attach(data)
y = natecon5
N = nrow(data)
Nstate = length(unique(data$stateID)) 
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "stateID"=stateID, "chgUnemp1615"= chgUnemp1615, "avgUnemp16"=avgUnemp16, 
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Nstate" = Nstate, "Ncat" = Ncat)

params = c("b.state", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "cut", "sigma.state")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2016_M1.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed=0508)
mcmc16_step1Model1_M1 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:21){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc16_step1Model1_M1)[1:21][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc16_step1Model1_M1)[order(mcmcID)] = c("chgUnemp1615", "avgUnemp16",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc16_step1Model1_M1, "mcmc16_step1Model1_M1")
rm(output); rm(mcmc16_step1Model1_M1)


## Model 2 in Table 3 ##
attach(data)
y = natecon5
N = nrow(data)
Nstate = length(unique(data$stateID)) 
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "stateID"=stateID, "chgUnemp1615"= chgUnemp1615, "avgUnemp16" = avgUnemp16, "chgInc1615"=chgInc1615, "inc16"=inc16,
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Nstate" = Nstate, "Ncat" = Ncat)

params = c("b.state", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "b.22", "b.23", "cut", "sigma.state")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2016.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed = 0508)

mcmc16_step1Model1_M2 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:23){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc16_step1Model1_M2)[1:23][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc16_step1Model1_M2)[order(mcmcID)] = c("chgUnemp1615", "avgUnemp16","chgInc1615", "inc16",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc16_step1Model1_M2, "mcmc16_step1Model1_M2")
rm(output); rm(mcmc16_step1Model1_M2)


## Model 3 in Table 3 ##
attach(data)
y = natecon5
N = nrow(data)
Ncounty = length(unique(data$countyID)) 
countyID = as.integer(droplevels(factor(countyID)))
Ncat = 5 # Categories in natecon5

data.jags = list("y"=y, "countyID"=countyID, "chgUnemp1615"= chgUnemp1615, "avgUnemp16" = avgUnemp16, "chgInc1615"=chgInc1615, "inc16"=inc16,
                 "younger30"=younger30, "older65"=older65, "female"=female, "black"=black, "latino"=latino, "noHighSchool"=noHighSchool, 
                 "HighSchool"=HighSchool, "someCollege"=someCollege, "fourCollege"=fourCollege, "income"=income, "unemployed"=unemployed, 
                 "fullTime"=fullTime, "partTime"=partTime, "ownHome"=ownHome, "newsInt"=newsInt, "ideol"=ideol,
                 "Dem"=Dem, "Rep" = Rep, "Ind"=Ind,
                 "N" = N, "Ncounty" = Ncounty, "Ncat" = Ncat)

params = c("b.county", "b.1", "b.2", "b.3", "b.4", "b.5", "b.6","b.7","b.8","b.9",
           "b.10", "b.11", "b.12", "b.13", "b.14", "b.15", "b.16","b.17","b.18","b.19",
           "b.20", "b.21", "b.22", "b.23", "cut", "sigma.county")
output = jags(data.jags, inits=NULL, params, model.file="Model1_2016_M3.bug", n.chains=1, n.burnin=20000, n.iter=70000, jags.seed = 0508)

mcmc16_step1Model1_M3 = as.mcmc(output)[[1]]
mcmcID=NA
for(i in 1:23){ # This is to adjust the order of the colnames of mcmc; the order is NOT 1,2,3,... BUT 1,10,100,101, ...
  ID=str_extract(colnames(mcmc16_step1Model1_M3)[1:23][i], "[0-9]+$") # Extract the digits located at the end of each string 
  mcmcID[i]=as.numeric(ID) 
}
colnames(mcmc16_step1Model1_M3)[order(mcmcID)] = c("chgUnemp1615", "avgUnemp16","chgInc1615", "inc16",
                                                   "younger30", "older65", "female", "black", "latino", "noHighSchool", 
                                                   "HighSchool", "someCollege", "fourCollege", "income", "unemployed", 
                                                   "fullTime", "partTime", "ownHome", "newsInt", "ideol",
                                                   "Dem", "Rep", "Ind") 
write.table(mcmc16_step1Model1_M3, "mcmc16_step1Model1_M3")
rm(output); rm(mcmc16_step1Model1_M3)


#############
# Table 6
#############

## Model 1 in Table 6 ##
model16 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp1615 + avgUnemp16 + 
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="stateID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc16_step1Model2_M1 = model16$mcmc[, 1:27]
colnames(mcmc16_step1Model2_M1)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp1615", "avgUnemp16",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc16_step1Model2_M1, "mcmc16_step1Model2_M1")

## Model 2 in Table 6 ##
model16 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp1615 +  avgUnemp16 + chgInc1615 + inc16 +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="stateID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc16_step1Model2_M2 = model16$mcmc[, 1:29]
colnames(mcmc16_step1Model2_M2)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp1615", "avgUnemp16", "chgInc1615", "inc16",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc16_step1Model2_M2, "mcmc16_step1Model2_M2")
rm(mcmc16_step1Model2_M2)


## Model 3 in Table 6 ##
model16 = MCMChlogit(fixed=pvote2 ~ VGood + Good + Bad + VBad +
                       chgUnemp1615 +  avgUnemp16 + chgInc1615 + inc16 +
                       younger30 + older65 + female + black + latino + noHighSchool + HighSchool + someCollege + fourCollege +  
                       income + unemployed + fullTime + partTime + ownHome +
                       Very + Few + ideol + Dem + Rep + Ind,
                     random=~1, group="countyID",
                     data=data, burnin=20000, mcmc=50000, thin=4,verbose=1,
                     seed=0508, beta.start=0, sigma2.start=1,
                     Vb.start=1, mubeta=0, Vbeta=1.0E6,
                     r=1, R=diag(1), nu=0.001, delta=0.001, FixOD=1)

mcmc16_step1Model2_M3 = model16$mcmc[, 1:29]
colnames(mcmc16_step1Model2_M3)=c("intercept", "VGood", "Good", "Bad", "VBad", 
                                  "chgUnemp1615", "avgUnemp16", "chgInc1615", "inc16",
                                  "younger30", "older65", "female", "black", "latino", "noHighSchool", "HighSchool", "someCollege", "fourCollege", "income", 
                                  "unemployed", "fullTime", "partTime", "ownHome",
                                  "Very", "Few", "ideol",  "Dem", "Rep", "Ind")
write.table(mcmc16_step1Model2_M3, "mcmc16_step1Model2_M3")
rm(mcmc16_step1Model2_M2)



###########################
#
# REPORT RESULTS
#
###########################

rm(list = ls())

quantileFun = function(x){quantile(x, c(0.025, 0.5, 0.975))}


## Table 1
mcmc08_step1Model1_M1 = read.table("mcmc08_step1Model1_M1")
mcmc08_step1Model1_M2 = read.table("mcmc08_step1Model1_M2")
mcmc08_step1Model1_M3 = read.table("mcmc08_step1Model1_M3")
mcmc08_step1Model1_M4 = read.table("mcmc08_step1Model1_M4")

quantile08_table1_M1 = apply(mcmc08_step1Model1_M1, 2, FUN=quantileFun)
quantile08_table1_M2 = apply(mcmc08_step1Model1_M2, 2, FUN=quantileFun)
quantile08_table1_M3 = apply(mcmc08_step1Model1_M3, 2, FUN=quantileFun)
quantile08_table1_M4 = apply(mcmc08_step1Model1_M4, 2, FUN=quantileFun)

round(quantile08_table1_M1[,c("chgUnemp0807", "avgUnemp08")], 3)
round(quantile08_table1_M2[,c("chgUnemp0807", "avgUnemp08","chgInc0807", "inc08")], 3)
round(quantile08_table1_M3[,c("chgUnemp0807", "avgUnemp08","chgInc0807", "inc08",
                              "avgGas", "totForcAugOct")], 3)
round(quantile08_table1_M4[,c("chgUnemp0807", "avgUnemp08","chgInc0807", "inc08",
                              "avgGas", "totForcAugOct")], 3)


## Table 4
mcmc08_step1Model2_M1 = read.table("mcmc08_step1Model2_M1")
mcmc08_step1Model2_M2 = read.table("mcmc08_step1Model2_M2")
mcmc08_step1Model2_M3 = read.table("mcmc08_step1Model2_M3")
mcmc08_step1Model2_M4 = read.table("mcmc08_step1Model2_M4")
mcmc08_step1Model2_M5 = read.table("mcmc08_step1Model2_M5")

quantile08_table4_M1 = apply(mcmc08_step1Model2_M1, 2, FUN=quantileFun)
quantile08_table4_M2 = apply(mcmc08_step1Model2_M2, 2, FUN=quantileFun)
quantile08_table4_M3 = apply(mcmc08_step1Model2_M3, 2, FUN=quantileFun)
quantile08_table4_M4 = apply(mcmc08_step1Model2_M4, 2, FUN=quantileFun)
quantile08_table4_M5 = apply(mcmc08_step1Model2_M5, 2, FUN=quantileFun)

round(quantile08_table4_M1[,2:7], 2)
round(quantile08_table4_M2[,2:9], 2)
round(quantile08_table4_M3[,2:11], 2)
round(quantile08_table4_M4[,2:11], 2)
round(quantile08_table4_M5[,2:10], 2)


## Table 2
mcmc12_step1Model1_M1 = read.table("mcmc12_step1Model1_M1")
mcmc12_step1Model1_M2 = read.table("mcmc12_step1Model1_M2")
mcmc12_step1Model1_M3 = read.table("mcmc12_step1Model1_M3")
mcmc12_step1Model1_M4 = read.table("mcmc12_step1Model1_M4")

quantile12_table2_M1 = apply(mcmc12_step1Model1_M1, 2, FUN=quantileFun)
quantile12_table2_M2 = apply(mcmc12_step1Model1_M2, 2, FUN=quantileFun)
quantile12_table2_M3 = apply(mcmc12_step1Model1_M3, 2, FUN=quantileFun)
quantile12_table2_M4 = apply(mcmc12_step1Model1_M4, 2, FUN=quantileFun)

round(quantile12_table2_M1[,c("chgUnemp1211", "avgUnemp12")], 3)
round(quantile12_table2_M2[,c("chgUnemp1211", "avgUnemp12","chgInc1211", "inc12")], 3)
round(quantile12_table2_M3[,c("chgUnemp1211", "avgUnemp12","chgInc1211", "inc12",
                              "avgGas", "totForcAugOct")], 3)
round(quantile12_table2_M4[,c("chgUnemp1211", "avgUnemp12","chgInc1211", "inc12",
                              "avgGas", "totForcAugOct")], 3)


# Table 5
mcmc12_step1Model2_M1 = read.table("mcmc12_step1Model2_M1")
mcmc12_step1Model2_M2 = read.table("mcmc12_step1Model2_M2")
mcmc12_step1Model2_M3 = read.table("mcmc12_step1Model2_M3")
mcmc12_step1Model2_M4 = read.table("mcmc12_step1Model2_M4")

quantile12_table5_M1 = apply(mcmc12_step1Model2_M1, 2, FUN=quantileFun)
quantile12_table5_M2 = apply(mcmc12_step1Model2_M2, 2, FUN=quantileFun)
quantile12_table5_M3 = apply(mcmc12_step1Model2_M3, 2, FUN=quantileFun)
quantile12_table5_M4 = apply(mcmc12_step1Model2_M4, 2, FUN=quantileFun)

round(quantile12_table5_M1[,2:7], 2)
round(quantile12_table5_M2[,2:9], 2)
round(quantile12_table5_M3[,2:11], 2)
round(quantile12_table5_M4[,2:11], 2)



## Table 3
mcmc16_step1Model1_M1 = read.table("mcmc16_step1Model1_M1")
mcmc16_step1Model1_M2 = read.table("mcmc16_step1Model1_M2")
mcmc16_step1Model1_M3 = read.table("mcmc16_step1Model1_M3")

quantile16_table3_M1 = apply(mcmc16_step1Model1_M1, 2, FUN=quantileFun)
quantile16_table3_M2 = apply(mcmc16_step1Model1_M2, 2, FUN=quantileFun)
quantile16_table3_M3 = apply(mcmc16_step1Model1_M3, 2, FUN=quantileFun)

round(quantile16_table3_M1[,c("chgUnemp1615", "avgUnemp16")], 3)
round(quantile16_table3_M2[,c("chgUnemp1615", "avgUnemp16","chgInc1615", "inc16")], 3)
round(quantile16_table3_M3[,c("chgUnemp1615", "avgUnemp16","chgInc1615", "inc16")], 3)


 # Table 6
mcmc16_step1Model2_M1 = read.table("mcmc16_step1Model2_M1")
mcmc16_step1Model2_M2 = read.table("mcmc16_step1Model2_M2")
mcmc16_step1Model2_M3 = read.table("mcmc16_step1Model2_M3")

quantile16_table6_M1 = apply(mcmc16_step1Model2_M1, 2, FUN=quantileFun)
quantile16_table6_M2 = apply(mcmc16_step1Model2_M2, 2, FUN=quantileFun)
quantile16_table6_M3 = apply(mcmc16_step1Model2_M3, 2, FUN=quantileFun)

round(quantile16_table6_M1[,2:7], 3)
round(quantile16_table6_M2[,2:9], 3)
round(quantile16_table6_M3[,2:9], 3)
